/* Program to calculate PI */

#include <stdlib.h>             /* Standard Lib    ISOC  */
#include <stdarg.h>             /* Variable args   ISOC  */
#include <stdio.h>              /* I/O lib         ISOC  */
#include <gmp.h>                /* GNU GMP Library       */
#include <math.h>


int main()
{
//double a=1, b, t=0.25, x=1, y;
mpf_t  a_mp, b_mp, t_mp, x_mp, y_mp, temp_mp, pi_mp;
unsigned int its, i;





//b = 1/sqrt(2);

//printf("\n Entre com o numero de casas decimais de precisao: >");
//scanf("%ud", &its);

//seta a precisao do ponto flutuante
mpf_set_default_prec(1024 * 40000);

mpf_init_set_d(a_mp,  1.0);
mpf_init_set_d(t_mp,  0.25);
mpf_init_set_d(x_mp,  1.0);
mpf_init_set_d(y_mp,  0.0);
mpf_init_set_d(temp_mp,  0.0);
mpf_init_set_d(pi_mp,  0.0);


mpf_sqrt_ui (temp_mp, 2);
mpf_ui_div (temp_mp, 1, temp_mp);
mpf_init_set(b_mp, temp_mp);

for(i=0; i<22; i++)
{
	//y = a;
	mpf_set (y_mp, a_mp);
	
	//a = (a + b) / 2;
	mpf_add (a_mp, a_mp, b_mp);
	mpf_div_ui (a_mp, a_mp, 2);

	//b = sqrt(b * y);
	mpf_mul (temp_mp , b_mp , y_mp);
	mpf_sqrt (b_mp, temp_mp);


	//t = t - x*((y - a)*(y-a));
	mpf_sub (temp_mp, y_mp, a_mp);
	mpf_mul (temp_mp, temp_mp, temp_mp);
	mpf_mul (temp_mp, x_mp, temp_mp);
	mpf_sub (t_mp, t_mp, temp_mp);


	//x = 2*x;
	mpf_mul_ui (x_mp, x_mp, 2);
}

//printf("\nPI = %.10000f\n\n", ((a+b)*(a+b))/(4*t));
mpf_add (temp_mp, b_mp, a_mp);
mpf_mul (temp_mp, temp_mp, temp_mp);
mpf_mul_ui (pi_mp, t_mp, 4);
mpf_div (pi_mp, temp_mp, pi_mp);


printf("\n");
mpf_out_str(NULL, 10, 10000000, pi_mp);
printf("\n");
return(0);
}
